function [narrowsnr, widesnr] = SNRbarplot(static,yaw,pitch,Fs,stimulusFreq,nHarmonics,excludeTrial)
%---input---
%staic/yaw/pitch：SSVEP数据，第一个维度为session，第二个维度为stimulus，第三个维度为trial，第四个维度为channel，第五个维度为samplepoints
%Fs：采样率
%stimulusFreq：刺激频率，e.g.[6.67 12.0 10.0 8.57]
%nHarmonics：谐波次数，用来计算wideSNR
%excludeTrial: 移除的trial，每一行表示一个trial，格式为[stimulus condition session trial]
%---output---
%narrowSNR：窄带信噪比
%wideSNR：宽带信噪比
narrowsnr = zeros(size(static,1),3,size(static,2));
widesnr = zeros(size(static,1),3,size(static,2));
for condition = 1:3
    switch condition
        case 1
            SSVEPdata = static;
        case 2
            SSVEPdata = yaw;
        case 3
            SSVEPdata = pitch;
    end
    for session = 1:size(SSVEPdata,1)
        for stimulus = 1:size(SSVEPdata,2)
            avl_trial_index = 0;
            nsnr_trial = [];
            wsnr_trial = [];
            for trial = 1:size(SSVEPdata,3)
                %移除错误数据
                avl_trial_flag = 1;
                for i = 1:size(excludeTrial,1)
                    excludei = excludeTrial(i,:);
                    if(all([stimulus condition session trial] == excludei))
                        avl_trial_flag = 0;
                        break;
                    end
                end
                if(~avl_trial_flag)
                    continue;
                end
                avl_trial_index = avl_trial_index + 1;
                epoch = SSVEPdata(session,stimulus,trial,:,:);
                epoch = reshape(epoch,size(SSVEPdata,4),size(SSVEPdata,5));
                [fftdata,~] = fast_fourier_transform(epoch',Fs);
                nsnr_trial(avl_trial_index) = narrow_signal_noise_ratio(fftdata, stimulusFreq(stimulus),Fs);
                wsnr_trial(avl_trial_index) = wide_signal_noise_ratio(fftdata, stimulusFreq(stimulus),nHarmonics,Fs);
            end
            narrowsnr(session,condition,stimulus) = mean(nsnr_trial);
            widesnr(session,condition,stimulus) = mean(wsnr_trial);
            display(['NarrowSNR: ',num2str(narrowsnr(session,condition,stimulus)), ' Session:', num2str(session),...
                ' Condition: ', num2str(condition), ' Stimulus: ', num2str(stimulus)]);
        end
    end
end

narrowsnrMean = reshape(mean(narrowsnr, 1),3,4);
widesnrMean = reshape(mean(widesnr, 1),3,4);
narrowsnrVar = reshape(var(narrowsnr, 0, 1),3,4);
widesnrVar = reshape(var(widesnr,0,1),3,4);
narrowsnrSEM = sqrt(narrowsnrVar / (size(narrowsnr,1)));
widesnrSEM = sqrt(widesnrVar / (size(widesnr,1)));
X = categorical({'6.67','8.57','10.0','12.0'});
X = reordercats(X,{'6.67','8.57','10.0','12.0'});
bar_plot(narrowsnrMean,narrowsnrSEM,X,'Narrow SNR of SSVEPs in different moving status and stimuli frequency',...
    'Stimulus Frequency [Hz]','SNR [dB]');
bar_plot(widesnrMean,widesnrSEM,X,'Wide SNR of SSVEPs in different moving status and stimuli frequency',...
    'Stimulus Frequency [Hz]','SNR [dB]');
set(gca,'YDir','reverse')
end



